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The direction dependence of the event rate in WIMP direct detection experiments provides a 
powerful tool for distinguishing WIMP events from potential backgrounds. We use a variety of (non- 
parametric) statistical tests to examine the number of events required to distinguish a WIMP signal 
from an isotropic background when the uncertainty in the reconstruction of the nuclear recoil direc- 
tion is included in the calculation of the expected signal. We consider a range of models for the Milky 
Way halo, and also study rotational symmetry tests aimed at detecting non-sphericity /isotropy of 
the Milky Way halo. Finally we examine ways of detecting tidal streams of WIMPs. We find that if 
the senses of the recoils are known then of order ten events will be sufficient to distinguish a WIMP 
signal from an isotropic background for all of the halo models considered, with the uncertainties 
in reconstructing the recoil direction only mildly increasing the required number of events. If the 
senses of the recoils are not known the number of events required is an order of magnitude larger, 
with a large variation between halo models, and the recoil resolution is now an important factor. 
The rotational symmetry tests require of order a thousand events to distinguish between spherical 
and significantly triaxial halos, however a deviation of the peak recoil direction from the direction of 
the solar motion due to a tidal stream could be detected with of order a hundred events, regardless 
of whether the sense of the recoils is known. 
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I. INTRODUCTION 

Weakly Interacting Massive Particle (WIMP) di- 
rect detection experiments aim to directly detect non- 
baryonic dark matter via the elastic scattering of WIMPs 
on detector nuclei, and are presently reaching the sen- 
sitivity required to detect neutralinos (the lightest su- 
persymmetric particle and an excellent WIMP candi- 
date). Since the expected event rates are very small ( 
0(10^^ — 1) counts kg~^day~^) distinguishing a putative 
Weakly Interacting Massive Particle (WIMP) signal from 
backgrounds due to, for instance, neutrons from cosmic- 
ray induced muons or natural radioactivity, is crucial. 
The Earth's motion about the Sun provides two potential 
WIMP smoking guns: i) an annual modulation 1] and 
ii) a strong direction dependence of the event rate. 
In principle the dependence of the differential event rate 
on the atomic mass of the detector (see e.g. Refs. |3,0|) 
is a third possibility, however this would require good 
control of systematics for detectors composed of different 
materials. 

The DAMA (DArk MAtter) collaboration have, us- 
ing ^ 100 kg of Nal crystals and an exposure of ~ 
1.1 X 10^ kgday, observed an annual modulation in their 
event rate Q, which they interpret as a WIMP signal. 
Taken at face value the allowed region of WIMP mass 
and cross-section parameter space is incompatible with 
the exclusion limits from other experiments (such as the 
Cryogenic Dark Matter Search, Edelweiss and Zeplin 
I) [g, however the form of the WIMP annual modula- 



tion signal, in particular its phase, depends sensitively 
on the local WIMP velocity distribution 0,11, H [l^ ^• 

It is, however, difficult to exclude the possibility that 
the annual modulation could be due to some previously 
unidentified background. 

The direction dependence of the WIMP scattering 
rate has several advantages over the annual modulation. 
Firstly, the amplitude of the annual modulation is typ- 
ically of order a few per-cent P, 3 while the event rate 
in the forward direction is roughly an order of magni- 
tude larger than that in the backward direction 0,1^. 
Secondly it is difhcult for the directional signal to be 
mimicked by backgrounds; in most cases (a point source 
in the lab is a possible exception) a background which is 
anisotropic in the laboratory frame will be isotropic in 
the Galactic rest frame as the, time dependent, conver- 
sion between the lab and Galactic co-ordinate frames will 
wash out any lab specific features. Furthermore the mean 
direction, in the lab, of WIMP induced recoils will vary 
over a sidereal day due to the rotation of the Earth . 
Designing a detector capable of measuring the directions 
of sub-100 keV nuclear recoils is however a considerable 
challenge. Prototype directional detectors based on roton 
anisotropy in liquid He 14], phonon anisotropy in BaF 



Other theoretical possibilities which could pe rha ps resolve this 
conflict in clud e spin-dependent interactions llll and inelastic 
scattering |12| . 
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crystals [l5l an d the anisotropic scintillation properties 
of stilbene | 16| | have been tested, however low pressure 
gas time projection chambers (TPCs), such as DRIFT 
(Directional Recoil Identification From Tracks) M,M 
and NEWAGE seem to offer the best prospects for 
a workable detector. In discriminating a WIMP signal 
from isotropic backgrounds it is therefore crucial to take 
into account the accuracy with which recoil directions 
can be measured by such a detector. Furthermore, the 
direction dependence of the recoil rate depends on the 
local WIMP velocity distribution |23| • Whilst this opens 
up the exciting possibility of experimentally probing the 
local dark matter distribution if WIMPs are detected (i.e. 
doing WIMP astronomy), it is also crucial that uncertain- 
ties in this distribution are accounted for when designing 
tests for discriminating a WIMP signal from isotropic 
backgrounds. 

Earlier studies have found that as few as 30 events 
would be required to distinguish a WIMP induced sig- 
nal from isotropic backgrounds [23, 0| and that with of 
order hundreds of events it may be possible to distin- 
guish between different models for the Galactic halo . 
In this paper we improve on these works by including 
the uncertainty in the reconstruction of the nuclear re- 
coil direction and applying non-parametric tests to un- 
binned data. In Sec. ^1 we describe our calculation of the 
nuclear recoil spectrum, including the modeling of the 
Milky Way halo fSec. lII A|l and the reconstruction of the 
recoil direction fSec. IIIlI|l . In Sec. IIIII we then apply an 
array of statistical tests aimed at probing the isotropy 

(Sec. nira 



rotational symmetry (Sec. IIII B|l and mean 
direction fSec lIIlC)) of a putative WIMP directional sig- 
nal, before concluding with discussion of our results in 
Sec. lIVI In Appendices A, B and C we outline the calcu- 
lation of the Earth's orbital velocity, the statistics used 
and the hypothesis testing formalism respectively. 



II. CALCULATING THE NUCLEAR RECOIL 
SPECTRUM 

A. Modeling the Milky Way halo 

The simplest possible model of the Milky Way (MW) 
halo is an isotropic sphere with density distribution p oc 
r~^, in which case the velocity distribution at all posi- 



tions within the halo is maxwellian: 
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where cto, the velocity dispersion, is related to the asymp- 
totic circular velocity (which we take throughout to be 
Vc = 220 kms"^) by ctq = ^/3/2vc. 

Observations and numerical simulations indicate that 
galaxy halos are in fact triaxial and anisotropic and con- 
tain substructure. Numerical simulations produce both 
prolate (c/6 > b/a) and oblate {c/b < b/a) halos (where 
(a, 5, c) are the lengths of the long, intermediate and 



short principal axes of the density distribution respec- 
tively) [23, H^. To date only cluster sized halos have 
been simulated in large numbers. Ref. jl^l found that 
the shape parameters of such large halos lie in the ranges 
c/a ~ 0.4 — 0.6 and b/a ~ 0.7 — 0.8. However the axial 
ratios vary with radius within a given halo, dependin g at 
least partly on the merger history of the halo. Ref. |2J| 
studied two high resolution simulations of MW size ha- 
los, one halo had c/a ~ b/a ^ 0.4, the other c/a ^ 0.5 
and b/a ~ 0.8. These simulations contain only collision- 
less matter, and the addition of baryons leads to halos 
that are less triaxial, especially in the central regions (251] . 
These ranges of axial ratios are in rough agreement with 
observational measurements of the shape of dark matter 
halos (for reviews see e.g. Ref. H^). 

The two MW-like halos in Ref. Q have velocity 
anisotropy at the solar radius (r = i?o ~ 8kpc) /3(i?o) ^ 
0.3 ±0.1 where (3 is defined as: 



P{r) = 1 - 



2(«2) 
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where (wg), (v^) and {v^) are the mean square velocity 
components in Galactic co-ordinates. 



1. Triaxial and anisotropic models 

We consider two anisotropic and/or triaxial halo mod- 
els: the logarithmic ellipsoidal model [23, and the 
Osipkov-Merritt model |29l l30l |. We summarize these 
models briefly below. For further details see either the 
original papers [Ulllls^l or Refs. |^and |33 |. 

The logarithmic ellipsoidal model|2^ |28| is the sim- 
plest triaxial, scale free generalization of the isothermal 
sphere, its shape and velocity anisotropy being indepen- 
dent of radius. This model was studied in detail, in the 
context of WIMP direct detection, by Evans, CaroUo and 
de Zeeuw '28^ and has since been used widely in calcula- 
tions of the expected WIMP signals ,8, _10., ^2] . In this 
model, it is assumed that the principle axes of the veloc- 
ity distribution arc aligned with conical co-ordinates, and 
in any of the planes of the halo conical co-ordinates coin- 
cide with local cylindrical polar co-ordinates and the lo- 
cal WIMP distribution can be approximated by a multi- 
variate Gaussian: 
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where {uB,,a^,az) are the radial, azimuthal and polar 
velocity dispersions respectively. These dispersions are 
given in terms of parameters p and q that are related to 
the halo axis ratios and 7 which is related to the isotropy 
parameter 3 |2ll33|. 

Osipkov-Merritt models |2^ Is^l provide self-consistent 
radially anisotropic velocity distribution functions for ha- 
los with spherically symmetric density profiles, with ve- 
locity anisotropy /3(r) = /(r^ + ''a) where is the 
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anisotropy radius. We follow Ref. |31l and assume that 
the MW has an NFW profile psj with scale radius 
20kpc ^ and fix Ta = 20 and 12kpc as in Ref. j cor- 
responding to P{Ra) = 0.14 and 0.31, and calculate the 
local velocity distribution using the fitting functions in 
Ref. 01. 

The parameters of the fiducial halo models which we 
consider are chosen to span the range of halo properties 
discussed above and are summarized in Table ^ Some 
of the triaxial models (numbers 6-9) are rather extreme 
(with (Tz << tJr, (T^), however they serve as a best /worse 
case scenario allow us to assess whether the deviation of 
the recoil spectrum from that produced by the standard 
halo is detectable/problematic. 



2. Tidal streams 

The extent to which the local WIMP distribution is 
smooth is an open question HE HE HZli however 
it is certainly possible that a stream of high velocity 
particles from a late accreted sub-halo could be pass- 
ing through the solar neighborhood 35, 36J and would 
produce a potentially distinctive signal in direct detec- 
tion experiments . Indeed one of the tidal tails of the 
Sagittarius dwarf galaxy (Sgr) passes close to the solar 
neighborhood jsEHEl- Helmi et al. jij previously dis- 
covered a group of halo stars in the solar neighborhood 
moving coherently with mean velocity in Galactic coor- 
dinatesujtr ~ (—65.0, 135.0, — 250.0) km s~^ and Freese 
et al. [421 have argued that these stars are in fact part 
of the Sgr tidal stream. This interpretation is somewhat 
controversial. Firstly (as noted by Freese et al.) there is a 
discrepancy between the measured velocity of these stars 
and that expected for the Sgr stream. Secondly, there is a 
second lower density stellar stream in the solar neighbor- 
hood moving in the opposite direction [4]| which is not 
expected for debris associated with Sgr. Finally the com- 
parison of detailed simulations of the kinematics of the 
Sgr stream with measured radial velocities of Sgr stream 
stars appears to indicate that the stream does not pass 
through the solar neighborhood (4^] . 

None the less there is a debris stream passing through 
the solar neighborhood regardless of whether or not it 
originated from the Sgr dwarf galaxy. We therefore fol- 
low Ref. [sE 113 s-nd assume that there is dark matter 
associated with the stream and model its velocity distri- 
bution of WIMPs as a maxwellian distribution with bulk 
velocity Vstr, and velocity dispersion astr- 

= (2./3)3/^.L [-^^) • 



^ The exact ch oice of profile has httle influence on the local velocity 
distribution |8l| : at the solar radius the slope is close to —2 
whatever profile is used. 
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TABLE I: Summary of the parameters of the fiducial halo 
models. 'SHM' denotes the standard halo model, 'LGE' the 
logarithmic ellipsoidal model with the Sun located on either 
the long (L) or intermediate (I) axis, 'OM' the Osipkov- 
Merritt model and 'SHM-I- Str' the standard halo model plus 
tidal stream contributing pstr/po of the local WIMP density, 
with bulk velocity and velocity dispersion as described in the 
text. In each case we take Vc = 220kms~^ and /3 is the 
anisotropy parameter, defined as in eq. 

We take a^u- = 30kms"\ and po = 0.07GeVcm"^(w 
0.25po)i (towards the lower and u pper ends of the range 
of values suggested in Ref. [sE 113 respectively and 
therefore providing a best case scenario for detectabil- 
ity), and for simplicity use the standard maxwellian ve- 
locity distribution [eq. (Q)] for the smooth background 
component. 



B. Modeling the detector response 

TPC based detectors are designed to measure the di- 
rections of nuclear recoils by drifting the ionisation pro- 
duced by recoils in the gas volume to a suitable charge 
readout plane. The difference in ionisation track lengths 
between electrons, alpha particles and nuclear recoils at 
low gas pressure allows efficient rejection of these back- 
grounds 0. To reduce charge diffusion along the track 
(which would restrict the spatial resolution and hence 
direction measurement and background discrimination) 
to thermal levels the detector can be filled with an elec- 
tronegative gas to give negative ion drift. The DRIFT 
collaboration have adopted CS2 as a suitable electroneg- 
ative gas; for further details see Ref. [TtIIT^ . 

The recoiling nucleus undergoes a large number of scat- 
terings (on average 30 (70) for a primary recoil energy 
of 20 (100) keV, with a gaussian distribution). These 
multiple scatterings together with the small, but finite, 
diffusion of drifted ionisation limits the accuracy with 
which the direction of the primary recoil can be recon- 
structed. In order to assess the likely accuracy of the 
track reconstruction we consider a TPC detector filled 
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FIG. 1: Simulated pixel read-out in the xy plane from a 40keV 
recoil. The strength of the signal in each pixel is obtained by 
integrating the voltage pulse for that pixel. The actual recoil 
track is shown by the solid line. The ellipse results from 
the moment analysis and its long axis gives the reconstructed 
direction. 



with 0.05 bar CS2 with a 200 /xm pitch micropixel readout 
plane, a 10cm drift length over which a uniform electric 
field of IkVcm"^ is applied. This gas pressure, mix- 
ture and electric field are chosen to match the design 
of the DRIFT-I detector 0. The micropixel readout 
is based on those described in Ref. 0, and the drift 
length is such that the rms diffusion over the full drift 
length is approximately equal to the pixel pitch "Zs"] . We 
use the SRIM2003 package |4fl] ^ to generate sulfur ioni- 
sation recoil tracks We simulate the drift and diffusion 
of the ionisation to the readout plane under the electric 
field and the subsequent generation of charge avalanches. 
The quenching factor (fraction of recoil energy going into 
ionization) is roughly 0.4 for a primary recoil energy of 
20 keV, in agreement with experimental data how- 
ever the direction reconstruction depends on the electron 
distribution and not their number. 

The charge distributions are then projected into the xy, 
yz and xz planes (as illustrated in Fig. The result- 
ing pixel hit pattern is approximately elliptical with long 
axis close to the primary recoil direction, which allows 
us to reconstruct the nuclear recoil directions via a mo- 
ment analysis. The distribution of the difference between 



^ This package was not designed to model recoils in gaseous tar- 
gets, however it accurately reproduces the recoil ranges in S 
found experimentally. 

* Carbon recoils are ignored as they account for less than 1 % of 
the recoil rate due to the dependence of the recoil rate and 
the low Carbon mass fraction. 



the primary recoil direction and the reconstructed track 
direction peaks at ~ 15°, decreasing weakly with increas- 
ing energy, and has a long large angle tail (the mean and 
root mean square deviation are ~ 25° and ^ 15°, again 
weakly decreasing with increasing energy). We include 
this stochastic uncertainty in our Monte Carlo simula- 
tions. 

The primary recoil energy threshold is 20 keV as below 
this energy the short track length ( 3-4 pixels) and mul- 
tiple scatterings make it impossible to reconstruct the 
track direction. With SRIM2003 generated recoils, we 
find near uniform distributions of ionisation along the 
tracks and therefore cannot determine the absolute signs 
of the reconstructed recoil vectors (i.e +x or — x). How- 
ever, we note that experimental measurements are really 
required to determine whether this absolute sign can be 
measured or not. We consider both possibilities below. 

We assume that the time-dependent conversion of re- 
coil directions measured in the laboratory frame to the 
galactic frame introduces no further errors in the mea- 
sured recoil directions due to, for instance, inaccuracy in 
the measured time of the event. 

We calculate the directional recoil rate above the 
20keV detector energy threshold via Monte Carlo sim- 
ulation so as to allow the inclusion of the detector res- 
olution effects discussed above. We work in Galactic 
{R, 0, z) co-ordinates, where R is directed towards the 
Galactic center, is in the direction of motion of the 
Local Standard of Rest (LSR) and z points towards 
the north Galactic pole. These co-ordinates are related 
to Galactic longitude I and latitude b by {R, </>, z) = 
(cos I cos b, sin / cos 6, sin 6) . 

Firstly the WIMP velocity distribution has to be trans- 
formed to the rest frame of the detector [/o(w^) 



/ffi('"'x)] Galilean transformation: 



The velocity of the Earth, w®(t), is the sum of the ve- 
locity of the LSR, Qlsr, the pecuhar velocity of the 
Sun, Vqp, and the Earth's orbital velocity about the 
Sun, Vorb{t). We take Qlsr = 220 kms"^ (see e.g. ^48]), 
— (10.0, 7.3, 5.2) kms~^ ^^'^ the results of 
the calculation of the Earth's orbital velocity detailed 
in Appendix A. The resulting WIMP velocity distri- 
bution in the detector rest frame is used to generate 
random incident WIMP velocities from the WIMP flux 
Uxf@{ux)i which then undergo isotropic scattering in the 
center of mass frame. We weight each recoil by a factor 
equal to the detector form factor (which is taken to have 
the Helm form appropriate for spin independent WIMP- 
sulfur scattering 4]) in order to take into account the 
variation of the scattering cross section with the momen- 
tum transfer to the nucleus. Finally we include the effects 
of multiple scattering on the reconstructed recoil direc- 
tion. 

The resultant, time averaged, WIMP and recoil distri- 
butions are shown as Hammer-Aitoff projections of the 
flux/rate in Galactic coordinates in Figs. [20 for mod- 
els 1 (standard halo model), 3 (logarithmic ellipsoidal 
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model with p = 0.9, q = 0.8 and /3 = 0.4 ) and 12 (stan- 
dard halo model with a 25% contribution to the local 
WIMP density from a tidal stream with properties as de- 
scribed in Sec. Ill A'2l) . For illustrative purposes we take 
Pa = 0.3 GeVcm"'^, — 100 GeV and for the recoil flux 
take the WIMP-proton elastic scattering cross-section to 
be a^p = 10~^ pb. Note the different scale for model 
12. We see that the WIMP flux produced by the triaxial 
halo model is significantly flattened with respect to the 
Galactic plane, however due to the finite detector resolu- 
tion the difference between the recoil rates predicted by 
models 1 and 3 is far smaller than the difference in the 
WIMP fluxes. For model 12, with the tidal stream, the 
peak direction of the WIMP flux and recoil rate both dif- 
fer significantly from the direction of motion of the Sun, 
suggesting that WIMPs associated with the stream of 
halo stars in the solar neighborhood could be detectable 
in a directional detector if the WIMP density is suffi- 
ciently large. 



III. APPLYING STATISTICAL TESTS 

A WIMP search strategy with a directional detector 
can essentially be divided into three regimes. Firstly, 
there is the simple search phase which aims to detect an 
anomalous recoil signal above that expected from back- 
grounds. Secondly, following the discovery of an anoma- 
lous recoil signal, there would be the confirmation stage. 
At this point the experiment would collect more data 
and aim to confirm the recoil signal as Galactic in origin 
by searching for the expected anisotropy in the recoil di- 
rections. In statistical terms the question posed at this 
point is 'Is the distribution of observed recoil directions 
isotropic?', more specifically 'How many events are re- 
quired to reject the hypothesis of isotropy (at a given 
confidence level)?'. Thirdly, assuming anisotropy is de- 
tected, the experiment would then collect more data and 
try and determine the form of the recoil anisotropy and 
hence the local WIMP velocity distribution. The halo 
models considered in Sec. Ill Al lead us to consider two 
simple hypotheses to test: 1) Does the distribution of 
recoil directions show evidence of flattening (i.e. a non- 
spherical halo)? and 2) Does the distribution of recoil 
directions show evidence for a signal from a tidal stream? 

As discussed in Sec. 13 previous work aimed at de- 
termining the number of events required to detect the 
anisotropy in the recoil directions has indicated that 
5 — 30 events are required to give a 90-95% rejection 
of isotropy |2^i21]. However, these analyses suffer from 
several non-optimal features. Neither analysis took into 
account the angular resolution of the recoil reconstruc- 
tion or the possibility that the absolute signs of the recoil 
vectors (i.e. their 'sense' +x or —x) may not be measur- 
able. In Ref. [23|, a Kolmogorov-Smirnov goodness of fit 
test was applied to the binned distribution of recoil cos 7 
values, where 7 is the angle between the recoil direction 
and the direction of motion of the Sun. Given the small 



+90° 




FIG. 2: Time averaged WIMP flux (top) and S Recoil rate 
above 20 keV (bottom) in Galactic {I, h) co-ordinates for halo 
model 1 and po = 0.3 GeV cm"'', = 100 GeV and a^p — 
10"^ pb. 

number of events expected at the confirmation stage of 
an experiment, binning data can lead to a significant loss 
of information and should therefore be avoided. Ref. 
used an unbinned likelihood analysis to determine the 
number of events required for a 95% confidence detec- 
tion of isotropy in 95% of experiments, for a range of 
halo models, and also the number of events required to 
distinguish a triaxial halo from the standard maxwellian 
halo model. The draw-back of likelihood analysis is that 
it requires the likelihood function for the recoil distribu- 
tion generated by a given halo model. It can therefore 
only give the relative likelihood of specific models (and 
parameter values), and there is no guarantee that any 
of the halo models considered is a good approximation 
to the local WIMP velocity distribution (especially since 
non-spherical and/or anisotropic halo models involve as- 
sumptions, which may or not be valid). As such this 
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FIG. 3: As Fig. 121 for model 3 (fogarithmic ellipsoidal model 
with p = 0.9, q = 0.8 and /3 = 0.4 ). 



analysis could not be applied to real data. 

Recoil directions constitute vectors, or, if the senses 
are unknown, undirected lines or axes, and so can equiv- 
alently be represented as points on a sphere. This al- 
lows us to use statistical inference methods developed 
for the analysis of spherically distributed data (for a re- 
view of this extensive field see the standard texts such 
as 113, 1^)- We investigate a variety of non-parametric 
statistics designed to test the isotropy fSecs. HTl All , ro- 
tational symmetry (Sec. IIII B|) and median direction 
(Sec. IIII C|l . When quoting results for the numbers of 
events required for detection, we focus on the benchmark 
halo models discussed in Sec. Ill Al however these tests do 
not make any assumptions about the form of the recoil 
spectrum and are hence equally applicable to any local 
WIMP velocity distribution. Recoil energies are not in- 
corporated in these statistics, but the energy provides no 
extra information about the degree of anisotropy and/or 
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FIG. 4: As Fig. El for model 12 (standard halo model with 
a 25% contribution to the local WIMP density from a tidal 
stream with properties given in Sec. Ill A"^ . 

rotational symmetry in the data. In addition, any analy- 
sis including the recoil energy would require (halo model 
dependent) assumptions about the form of the recoil en- 
ergy spectrum as a function of direction. 

Throughout we assume zero background. This is a rea- 
sonable expectation for experiments such as DRIFT-II 
made from low activity materials with efficient gamma 
rejection and shielding [52]. However, non-zero back- 
grounds can be incorporated into the statistical tests. 
For real data, one would simply calculate the value of 
the appropriate test statistic using the observed sam- 
ple of recoil vectors/axes and use its null distribution 
(Appendix Bl) to determine the probability of the sam- 
ple being isotropic. No assumption about the level of 
background contribution to the sample is neccessary, but 
there is an implicit assumption that background event 
directions are isotropically distributed. 

In the case of our simulations, it would be neccesary 
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to fix the background rate and exposure so that the ex- 
pected number of background events Ni, is known. A 
variable number of WIMP events can then be added 
to the A'b background events, and the formahsm of Ap- 
pendix C used to determine the value of required to 
reject isotropy of the N^! + Nh recoil sample at a given 
confidence level. Since the background rate and expo- 
sure are dependent on the experiment, investigating the 
effect of varying background rates is beyond the scope of 
this paper. Therefore, and given that the next genera- 
tion of directional detectors are expected to have essen- 
tially zero background [52l |. we assume zero background 
to provide benchmark figures. Nevertheless, the formal- 
ism presented in this paper could be used to predict the 
directional sensitivity of a detector with non-zero back- 
ground once the background rate and exposure have been 
estimated. 



A. Tests of isotropy 

The first fundamental question posed by a directional 
detector observing an anomalous recoil signal is 'Are the 
observed recoil directions consistent with an isotropic dis- 
tribution?'. The most general tests of isotropy are those 
which do not depend on the coordinate system in which 
the sample vectors/axes are measured. This indepen- 
dence from the coordinate system means no assumption 
about the form of the anisotropy in the signal is required. 
While this is usually an advantage, coordinate indepen- 
dent tests cannot make use of any information that is 
available about the expected form of the anisotropy. If 
the local WIMP distribution is smooth then the year av- 
eraged recoil flux will be peaked in the direction of the 
Sun's motion at (;0,5o) = (87.5°, 1.3°). 

For each of the statistics discussed in Appendix IB II 
we determine the minimum number of events required to 
reject isotropy of recoil directions at 90 (95)% confidence 
in 90 (95)% of experiments (i.e. for rejection and accep- 
tance probabilities of Rc = Ac = 0.9 and 0.95), iViso, 
as described in Appendix lUl for each statistic and halo 
model considered, including the detector response in the 
calculation of the recoil distributions. For the standard 
maxwellian halo model we also carry out the calculation 
neglecting the reconstruction uncertainty and for a hypo- 
thetic perfect detail which also has zero energy threshold 
, in order to assess the effect of the reconstruction uncer- 
tainties and the finite energy threshold on the number of 
events required. The results are tabulated in Table HTl for 
each of the fiducial halo models for each statistic. 

Overall, the two coordinate dependent tests of isotropy, 
{cos 9) and (| cos6'|), have the lowest A^iso values and are 
thus the most powerful tests for rejecting isotropy. All 
the coordinate independent tests typically require 1.5- 
2 times more events to reject isotropy for a given halo 
model (for both vectorial and axial data). This is per- 
haps not surprising for the models 1-11 however, the co- 
ordinate dependent tests are also the most powerful for 



rejecting isotropy in model 12 which includes a 25% con- 
tribution to the local WIMP density from a tidal stream. 

The vectorial coordinate independent tests (W*, A and 
J-) all have very similar values of A'iso and are hence 
equally powerful for rejecting isotropy of recoils. The 
same is true of the axial coordinate independent tests 
{B* and Q), however A^iso is of order ten times larger for 
these tests than for the vectorial tests, and there is much 
greater variation between the fiducial halo models. Both 



Halo 




for (i?c,Ac) 


= (0.90, 0.90) 


Model 


Vectorial Statistics 


Axial Statistics 




W* 


A 




(cos 6) 


B* 


G 


(Icos^l) 


1 


12 


12 


13 


7 


165 


167 


81 


1 (no) 


10 


10 


10 


5 


83 


84 


40 


1 (per) 


18 


18 


18 


10 


531 


538 


258 


2 


12 


12 


12 


7 


114 


114 


57 


3 


14 


14 


15 


8 


157 


159 


93 


4 


12 


12 


13 


7 


149 


151 


74 


5 


14 


14 


15 


8 


157 


159 


93 


6 


11 


11 


11 


6 


67 


67 


36 


7 


14 


14 


14 


8 


88 


90 


57 


8 


13 


13 


14 


7 


175 


178 


86 


9 


15 


15 


16 


9 


264 


267 


146 


10 


15 


15 


16 


8 


280 


284 


149 


11 


12 


12 


12 


7 


125 


126 


62 


12 


14 


14 


14 


8 


221 


222 


159 






for (7?c,Ac) 


= (0.95,0.95) 


1 


18 


18 


19 


11 


235 


235 


131 


1 (no) 


13 


13 


14 


9 


120 


120 


65 


1 (per) 


25 


25 


26 


15 


767 


769 


426 


2 


17 


17 


18 


10 


161 


163 


93 


3 


20 


20 


21 


12 


225 


224 


152 


4 


18 


18 


19 


11 


215 


214 


120 


5 


20 


20 


21 


12 


222 


22 


153 


6 


16 


16 


16 


10 


96 


97 


59 


7 


19 


20 


20 


12 


125 


125 


93 


8 


18 


18 


19 


11 


252 


252 


142 


9 


21 


21 


22 


13 


378 


383 


241 


10 


21 


21 


21 


12 


402 


405 


247 


11 


17 


17 


17 


10 


180 


182 


102 


12 


19 


19 


20 


12 


316 


320 


264 



TABLE II: Number of recoil events required to reject isotropy 
of recoil directions, A''iso, at 90 (95)% confidence in 90 (95)% of 
experiments for each test statistic and halo model considered, 
including the uncertainty in the nuclear recoil reconstruction 
and for the standard halo model ignoring the uncertainty in 
the nuclear recoil reconstruction (denoted by 'no' in the ta- 
ble) and for a hypothetical perfect detector with zero thresh- 
old and no reconstruction uncertainities (denoted by 'per'). 
The test statistics have been divided into those suitable for 
vectorial data and those applicable to axial data (i.e. when 
the senses of the recoil directions are not known). 
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of these effects are due to the form of the recoil direction 
distribution and the way in which lack of knowledge of 
the recoil sense changes this form. While most of the 
recoil arrival directions lie in the forward hemisphere, 
< ^ < 180°, a non-negligible (halo model dependent) 
number lie in the backward hemisphere. For the vecto- 
rial tests, these backward events are not a problem as the 
tests can distinguish between the forward and backward 
hemispheres. When only the recoil axis is known, events 
in the backward hemisphere are effectively 'added' to the 
forward hemisphere, reducing the degree of anisotropy. 
Halo model 6 has the smallest rate in the backward hemi- 
sphere and thus the lowest TVjso for both vectorial and 
axial tests, while model 10 has the largest backward rate 
and therefore has higher Niso values with axial tests. In 
general the broader the local WIMP velocity distribution, 
the greater the number of events in the backward hemi- 
sphere, and hence a larger number of events are required 
to reject isotropy if the sense of the recoil directions are 
not known. 

Ignoring the detector response only marginally reduces 
the number of events required for the vectorial tests, how- 
ever the number of events required for the axial tests are 
decreased by a factor of roughly two. For the hypotheti- 
cal detector with perfect reconstruction and zero energy 
threshold the number of events required increases, rela- 
tive to the 20 keV perfect reconstruction case, by a factor 
of ^ 6) for the vectorial (axial) statistics. The num- 
ber of events required increases since, as is well known , 
the anisotropy of the recoil directions decreases with de- 
creasing recoil energy. The total recoil rate with a keV 
threshold is roughly twice that with a 20 keV threshold, 
therefore for the vectorial statistics the exposure required 
to detect a WIMP signal is essentially unchanged, while 
for the axial statistics (i.e. if it is not posible to mea- 
sure the senses of the recoils) the lower threshold would 
actually increases the exposure required. 

We now translate the number of events required to re- 
ject isotropy for the realistic detector into the equivalent 
detector exposure, E, required to observe this number of 
events. If the senses of the recoil directions are observed, 
isotropy could be rejected at 95% confidence in 95% of 
experiments for ao 5x 10~^ pb and po 0.3 GeVcm"^ 
with an exposure of ~ 10^ kg day (i.e. a 100 kg detec- 
tor operating for a period of 2-3 years). For a detector 
only capable of measuring the recoil axes a 10^ kg day ex- 
posure would be able to reject isotropy at the same con- 
fidence level and acceptance down to ctq ^ 3 x 10~^pb. 



B. Tests for rotational symmetry 

Once an observed anomalous recoil signal is found to 
be incompatible (at some confidence level) with isotropy 
the next question to pose is 'Is the observed distribu- 
tion consistent with a spherical halo?'. A generic feature 
of triaxial halo models is a flattening of the recoil dis- 
tribution towards the galactic plane, or more generally. 



flattened along one principal axis of the local velocity dis- 
tribution. This type of distribution can be probed using 
tests for rotational symmetry about a specified direction 
or axis. In the case of smooth halo models this direc- 
tion/axis is the direction of motion of the Earth through 
the MW halo, which over the year, averages to the direc- 
tion of motion of the Sun, (Zq, 6©) ^. 

We focus here on the two fiducial halo models for which 
the recoil distribution deviates most from that produced 
by the standard halo model: the logarithmic ellipsoidal 
models 5 {p = 0.9, q = 0.8, /3 = 0.4, long axis) and 
7 {p = 0.72, q = 0.7, /3 = 0.4, intermediate axis) 
and find the number of events required to reject rota- 
tional symmetry, A'^rot, using the Kuiper statistic V* as 
defined in Appendix IB 21 We find that halo model 5 has 
A < 0.4 over a large range of A^ values, and of order 
5000-8000 events would be required to reject rotational 
symmetry at 90 -f % confidence with similar acceptance. 
In contrast for halo model 7 (which is more extreme) 
N,ot = 1170(1710) for ^ ^ 0.9(0.95). So while 
rotational symmetry of the recoil distribution can be re- 
jected at high confidence and acceptance for this, rather 
extreme, halo model, it requires 10-50 times the num- 
ber of events required to reject isotropy (depending on 
whether vectors or axes are measured in the latter case). 

Converting these numbers to the equivalent exposures 
in kg day gives E{Rc = Ac = 0.9(0.95)) = 7.8(11) x 
IO^'^ctq^^Pq ^. Hence for the 10^ kg day exposure con- 
sidered in the previous subsection, rotational symme- 
try could be rejected at 90% confidence down to (Jq ^ 
3 X 10"'' pb for Po 0.3GeVcm~'^, independent of 
whether the senses of the recoils can be measured. The 
sensitivity of a TPC-based directional detector for re- 
jecting a spherical halo if the real MW halo is signif- 
icantly triaxial using this statistic, is therefore several 
orders of magnitude lower than for rejecting isotropy of 
recoil directions. Devising other tests that may be a more 
powerful probe of the flattening of the recoil distribution 
produced by non-spherical halo models is therefore an 
important task. 



C. Tests for mean direction 

Tidal streams will typically have a velocity disper- 
sion which is small compared to their velocity relative 
to the solar system and therefore the recoil distribu- 
tion due to WIMPs from a stream will be peaked in 
the hemisphere whose pole points in the direction of 
the stream velocity. The net (stream plus smooth back- 
groimd WIMP distribution) peak direction depends on 



^ For experiments with an exposure that is non-uniform in time, 
the mean direction could be calculated by averaging the Earth's 
velocity vector over the, time-dependent, exposure. 

^ More spherical models that produce less flattened recoil distribu- 
tions would require more events to reject rotational symmetry. 
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FIG. 5: Variation of the angle (in degrees) between the mean 
recoil directions of the standard halo and the standard halo 
plus a tidal stream with Vstr ~ (—65.0, 135.0, —250.0) kms~^ 
function of the fraction of the local WIMP density contributed 
by the stream, pstr/po- 



exposure would be sufficient to identify the presence of 
a tidal stream with properties as described in Sec. Ill A 21 
comprising 25% of the local WIMP density for WIMP- 
proton cross sections down to ctq ^ 8 x 10~® (1 x 10"^) pb 
and po — 0.3GeVcm~'^. 

We caution that these numbers have been calculated 
for parameter values at the optimistic ends of the ranges 
estimated in Ref. [sl, i-^- high density and low 
velocity dispersion. A lower stream density and/or a 
higher velocity dispersion would give a peak recoil direc- 
tion closer to the mean direction of motion of the Sun, 
and make the deviation due to the stream harder to de- 
tect. For the most pessimistic case of Pstr/po = 0.003 
and (Tstr — 70kms~^, the net WIMP flux and recoil rate 
are virtually indistinguishable from those for the stan- 
dard halo model and tens of thousands of events would 
be needed to detect the stream. In general the directional 
detectability of cold streams of WIMPs will clearly de- 
pend on how much their bulk velocity deviates from the 
direction of solar motion. 



the fraction of the local density that the tidal stream 
contributes. Fig. |31 shows the angle between the direc- 
tion of the Sun's motion (^0,60) — (87.5°, 1.3°) and 
the direction of the peak in the recoil distribution pro- 
duced by a maxwellian halo plus a stream with velocity 
Vsti- = (—65.0, 135.0, —250.0) kms""'^ as a function of the 
local WIMP density contributed by the stream. For the 
Pstr/£n ~ 1 — 30% range of stream densities suggested 
by [33, tbe angle is in the range 2 — 15°. This fea- 
ture suggests that the presence of a WIMP stream in the 
solar neighborhood could be tested for by comparing the 
median direction of the observed recoil vectors with the 
(known) direction of solar motion. 

For Pstr/po = 25% (which is at the upper end of the 
range of values suggested in [3^ j and therefore pro- 
vides a lower limit on the number of events required), 
we find that the minimum number of events required to 
reject (^0,60) as the median direction, A^dir, using the 
test as defined in Appendix El is 201 (294) for i?c = 
Ac = 0.9 (0.95). Comparing these numbers with those re- 
quired to reject isotropy of recoils for this stream density, 
it can be seen that ^ 10—20 times more events are needed 
to reject the direction of the solar motion as the me- 
dian recoil direction. These numbers of events would be 
achievable with a 10^ kg day exposure for WIMP-proton 
cross sections down to ctq ~ 4 x 10~® (6 x 10^^) pb for 
Po = 0.3GeVcm-^ 

In the case of axial recoil directions, the test cannot 
be applied, however, the rotational symmetry test pre- 
sented in Appendix lB 2l can be used since the recoil distri- 
bution will not be rotationally symmetric about {Iq, Bq) 
in the presence of a stream and we find A^rot = 390 (574) 
for i?c = Ac = 0.9(0.95). These numbers are roughly 
twice those for the test, and also for rejecting isotropy 
with axial data. Thus with only axial data, a 10^ kg day 



IV. DISCUSSION 

We have studied the application of non-parametric 
tests, developed for the analysis of spherical data [50ll5ll |. 
to the analysis of simulated data as expected from a 
TPC-based directional WIMP detector 0,0], taking 
into account the uncertainties in the reconstruction of 
the nuclear recoil directions. These tests, unlike likeli- 
hood analysis, do not require any assumptions about the 
form of the local WIMP velocity distribution. This is ad- 
vantageous as even if the properties (i.e. shape, velocity 
anisotropy, and density profile) of the MW halo at the so- 
lar radius were accurately determined there would likely 
be a wide range of local velocity distributions consistent 
with these properties. 

For a range of fiducial anisotropic and/or triaxial halo 
models, with parameters chosen to reproduce the range 
of properties found in simulations/observations of dark 
matter halos, we calculate the number of events re- 
quired to distinguish the WIMP directional signal from 
an isotropic background using a variety of tests. The 
most powerful test is the, co-ordinate system dependent, 
test of the mean angle between the observed nuclear re- 
coils and the direction of motion of the Sun, which takes 
advantage of the fact that, for a smooth halo, the re- 
coil rate is peaked in the direction of the Sun's motion. 
We find that if the senses of the recoils are known then 
of order ten events will be sufficient to distinguish a 
WIMP signal from an isotropic background for all of the 
halo models considered, with the uncertainties in recon- 
structing the recoil direction only mildly increasing the 
required number of events. If the senses of the recoils 
are not known the number of events required is an or- 
der of magnitude larger, with a large variation between 
halo models and the recoil resolution is now an impor- 
tant factor. The number of events required would be 
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significantly larger if the WIMP velocity distribution in 
the rest frame of the detector is close to isotropic |20f . 
which could be the case if the MW halo is co-rotating or 
if the local dark matter density has a significant contri- 
bution from a cold flow with direction so as to cancel out 
the front-rear asymmetry from the smooth background 
distribution [s^ , however for halos formed hierarchically 
(as is the case in Cold Dark Matter cosmologies) nei- 
ther of these possibilities are expected to occur (see e.g. 
Ref. [3 113 )• Our conclusions for the case where the 
recoil senses are known are in broad agreement with the 
results of Ref. . It is reassuring that the realistic mod- 
eling of detector properties and the use of statistics that 
can be applied to real data do not significantly degrade 
the expected detection power of directional experiments. 

Halo models which produce significantly different 
WIMP flux distributions, unfortunately give very simi- 
lar recoil rates. We find that distinguishing between halo 
models, in particular determining whether the MW halo 
is (close to) spherical, using tests of rotational symmetry 
will require thousands of events. 

Gondolo has shown that the recoil momentum spec- 
trum is the Radon transform of the WIMP velocity dis- 
tribution [s^ . In theory the transform could be inverted 
to directly measure the WIMP velocity distribution from 
the measured recoil momentum distribution, however the 
inversion algorithms available in the literature require 
large numbers of events |54j . Furthermore, as we have 
shown, the uncertainties in the reconstruction of the re- 
coil directions due to multiple scattering and diffusion 
have a significant effect on the observed recoil distribu- 
tion and it is not possible to apply this inversion to real 
data. Developing techniques to distinguish between the 
recoil distributions expected from different halo models, 
for finite amounts of data and taking into account the 
experimental resolution, is therefore of key importance 
for the full astronomical exploitation of data from direc- 
tional detectors. 

If a significant fraction of the local dark matter is in 
a cold flow (with velocity dispersion far smaller than its 
bulk velocity) then the peak recoil direction will deviate 
from the direction of the solar motion. In fact a stream of 
halos stars has been found in the solar neighborhood ^Jl 
and using values for the properties of the stream at the 
optimistic end of the ranges estimated in Ref. [H, 
this deviation could be detected with of order a hundred 
events even if the sense of the nuclear recoils is not mea- 
surable. This illustrates that if the Galactic dark matter 
is in the form of WIMPs then WIMP directional detec- 
tors, such as DRIFT [Hill, will be able to do 'WIMP 
astronomy'. 
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APPENDIX A: ORBITAL VELOCITY OF THE 
EARTH 

The Earth's orbit around the Sun lies in the Eclip- 
tic plane and is approximately circular with eccentric- 
ity e = 0.0167 and semi-major axis a = 1.000 AU= 
1.496 X 10^^ m. The tangential and radial components 
of the Earth's orbital velocity at a given time t are cal- 
culated from ISal 



vrit) 

VR{t) 

v{t) 



27ra(l + e cos v{t)) 
2T:aes\Tiv{t) 
Aq (t) — tu — tt , 



(Al) 



where XQ{t) is the ecliptic longitude of the Sun, m — 
103° is the argument of perihelion of the Earth and 
P is the orbital period. The Solar ecliptic longitude at 
time t, accurate to 0.01° for t between 1950 and 2050, 
is calculated from XQ{t) = 280.460 4- 0.9856474i -|- 
1.915 sin (7 -1-0. 020 sm2g where t is the time given in terms 
of the Juhan Day JD via t = JD - 2451545.0 and g is 
the mean anomaly of the orbit given hy g = 357.528 -I- 
0.9856003i. Using Eqs. EH the Earth's orbital velocity 
in the Ecliptic coordinate system is calculated as 



where 



Ait) 
B[t) 



A{t)vR{t) - B{t)vT{t) 
A{t)vT{t) ~ B{t)vR{t) 




cos ^{t) COS w — sin ^{t) sin w , 
sin v{t) cos vj — cos ^{t) sin w . 



(A2) 
(A3) 



(A4) 
(A5) 



This velocity is then transformed to the galactic coordi- 
nate system via the rotations 



i;orb(i) = GEi;^,b(t), 
with rotation matrices 
(l 

E = 



(A6) 



(A7) 



where e = 23.5° is the obliquity of the ecliptic plane to 
the equatorial plane of the Earth and [ ssf 

/ -0.054876 -0.873437 -0.483835 \ 
G = 0.494109 -0.444830 0.746982 . (A8) 
\ -0.867666 -0.198076 0.455984 / 




11 



APPENDIX B: STATISTICAL TESTS 
1. Tests of isotropy 

The simplest coordinate independent statistic for vec- 
torial data is the so-called modified Rayleigh- Watson 
statistic W* ^ . For a sample of N miit vectors Xi, W* 
is defined as IsSillEi 



1 



1 
2N 



1 



107V 



(Bl) 



where W is the (unmodified) Rayleigh- Watson statistic 

3 



and TZ is the Rayleigh statistic: 

N 



n = 



(B2) 



(B3) 



The value of W* becomes larger as the degree of 
anisotropy increases and for isotropically distributed vec- 



tors, yV* is asymptotically distributed as xi |57ll58| . We 
find, using Monte Carlo simulations, that the difference 
between xi and the true distribution of W* for isotropic 
vectors in the large W* tail of the distribution is less 
than 2% for N > 30. For iV < 30 the xl distribution 
significantly underestimates the true probability distri- 
bution and we calculate the probability distribution from 
the exact probability distribution of 7?., as described in 
Ref. 

The Rayleigh- Watson statistic has the drawback that 
it is not sensitive to distributions which are symmetric 
with respect to the center of the sphere and therefore 
can not be used with axial data. The Bingham statistic 
avoids this problem and is based on the scatter (or 
orientation) matrix T of the data, defined as [sol . fssL Isof 



1 ^ 
N ^ 



4=1 



XiXi Xiyi XiZ'i 

Vi'^i ViVi yi-^i 
Z'lX'l Z'l'tJ'l Z^Zl 



(B4) 



where (xi, yi, Zi) are the components of the i-th vector or 
axis. This matrix is real and symmetric with unit trace, 
so that that the sum of its eigenvalues (fc = 1, 2, 3) is 
unity, and for an isotropic distribution all three eigenval- 
ues should, modulo statistical fiuctuations, be equal to 
1/3. Bingham's modified statistic 5* jMElil 



= B l- 



1 

N 



47 

84 



13 

147" 



B 



5292 



:B' 



(B5) 



This modified statistic, and the others considered in this section, 
have the advantage of approaching their large N asymptotic dis- 
tribution for smaller A'^ than the unmodified statistic. 



where B is the (unmodified) Bingham statistic 



B 



15N 



E 

fe=i 



(B6) 



measures the deviation of the eigenvalues from the 
value expected for an isotropic distribution. For isotrop- 
ically distributed vectors/axes B* is asymptotically dis- 
tributed as and we have found, via Monte Carlo simu- 
lations, that the difference between the underlying prob- 
ability distribution and Xs in the large B* tail of the 
distribution is always less than 10% and is less than 2% 
for N > 10. 

The most general test of the uniformity of a sample 
of A'^ unit vectors or axes is provided by the statistics of 
Beran ^31 and Gine 163 , which are defined as 



A 



N - 



^ N-l N 

A^E E V'..' 

^ = T - /v^E E ^i^'^' 



(B7) 
(B8) 



where ipij is the angle between the i-th and j-th direc- 
tions. Beran's statistic A tests for distributions which are 
asymmetric with respect to the center of the sphere, and 
Gine's statistic G tests for distributions which are sym- 
metric with respect to the center of the sphere. A suitable 
statistic for testing uniformity against all possible alter- 
natives is therefore the combination = A + G |51] . In 
the case of axial data G alone can be used. The asymp- 
totic distributions of ^, G and have not, to date, been 
calculated. We therefore generate the probability distri- 
butions for these statistics under the null hypothesis via 
Monte Carlo simulation. 

A suitable coordinate system dependent statistic which 
uses the fact that, for smooth halo models, the WIMP 
recoil distribution is expected to be peaked about 



N 



for vectorial data. 



N 



(B9) 



(BIO) 



for axial data, where 6i is the angle between (^0,60) 
and the ith vector/axis. For N isotropic vec- 
tors (cos 0) ((I cos 0|)) can take values on the interval 
[—1,1] ([0, 1]) and, due to the central limit theorem, has 
a gaussian distribution with mean (0.5) and variance 
1/3 (l/3iV) [g^]. The larger the concentration of recoil 
directions towards (Z0 , Bq ) the larger these statistics will 
be. We calculate the probability distribution function 
of (cos 6*) and (|cos0|) as a function of N for the null 
hypothesis of an isotropic recoil spectrum using Monte 
Carlo simulations. 
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2. Tests for rotational symmetry 

A test of rotational symmetry about some hypothe- 
sized direction (^o, "^o) (vahd for vectors or axes) can be 
performed by first rotating the sample vectors or axes so 
that their polar angles are measured relative to (^Oi^o)- 
The resultant azimuthal angles are then divided by 27r 
{Xi = (j)i/2n) and sorted in ascending order. For rota- 
tional symmetry, the distribution of the ordered normal- 
ized angles should be uniform between and 1. This hy- 
pothesis can be tested by using the Kuiper statistic 51| , 
which is related to the well-known Kolmogorov-Smirnov 
test, but has the advantage of being invariant under cyclic 
transformations and equally sensitive to deviations at all 
values of X. The modified Kuiper statistic is defined as 



V* = V i N^/^ + 0.155 



0.24 



(Bll) 



where V is the (unmodified) Kuiper statistic 

V^V++V-, (B12) 

and 



I 

N 



maxi— -X, ), i = l,...,N (B13) 



max X 



i - 1 
N 



(B14) 



As there is no general formula for the distribution of V"^ 
under the null hypothesis of rotational symmetry, we use 
Monte Carlos to generate the null probability distribu- 
tion assuming the standard halo model ^. 



3. Tests for mean direction 

For vectorial data, the median direction is defined 
as the direction {9, <j)) which minimizes the sum of ar- 
clengths between {9, (f>) and the sample vectors [HJ, and 
can found by minimizing the quantity 



M 



N 



VVi 



ZZi) . 



(B15) 



while for axial data an estimate of the principal axis 
of the distribution is provided by the eigenvector corre- 
sponding to the largest eigenvalue of the scatter matrix 
T. A test of compatibility between the measured median 
direction and a hypothesized median direction {9o,(f)o) 
for vectorial data can be performed .51] by first rotating 



the data so that the polar angles are measured relative 
to a pole in the direction of the sample median {9, (f>) and 
then calculating the matrix S 



S = 



1 / CTll cri2 

2 \ (721 t722 



where 



1 V- 

an = 1 + — 2^cos2(/), 

1=1 

1 ^ 

CT22 = 1 - XI 2*^^ 



1 ^ 

CT12 ^21 = X! ™ 2"^* 



(B16) 

(B17) 
(B18) 

(B19) 



from the rotated sample vectors. The sample vectors are 
rotated again so that their polar angles are measured 
relative to a pole in the direction of the hypothesized 
median using and the vector 



C7 = Ar-i/2 ( Ecos(/.o 



(B20) 



where (/)^ is the azimuthal angle of the ith sample vector 
relative to a pole at (0O;0o) is calculated. Finally the 
test statistic X^, is defined as 



(B21) 



and is distributed as X2 f^^i' samples drawn from a distri- 
bution with median direction (^Qj 4>a)- 



APPENDIX C: HYPOTHESIS TESTING 

For each halo model considered we use Monte Carlo 
simulations to generate N recoil scattering events in each 
of 10^ experiments, for values of N between 5 and 400 
(the lower value corresponding to the point at which an 
anomalous recoil signal would first be identified at high 
confidence). For each experiment the test statistic T is 
calculated from the N recoil directions used to give the 
probability distribution of the statistic, pi{T;N), as a 
function of the number of recoil events. The probability 
distribution for the null hypothesis of isotropy, Pq{T; N), 
is calculated using analytical expression where available 
and otherwise via Monte Carlo simulations. As in stan- 
dard hypothesis testing, the overlap between these two 
distributions allows the probability ^ with which the null 
and alternative hypotheses can be rejected or accepted 



This may appear to introduce a model dependence into the test, 
however any other distribution rotationally symmetric about 

(/0, 60) would give the same null distribution for V* . ^ We use the frequentist definition of 'probability' throughout. 
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to be calculated. For a given value 7^, the rejection fac- 
tor R is the probability of measuring T < 7^ if the null 
hypothesis is true: 

R = [ 'po{T;N)dT. (CI) 
Jo 

The rejection factor thus gives the confidence level with 

which the null hypothesis can be rejected given a par- 
ticular value of the test statistic T = %. For the same 
value of the test statistic Tg, the acceptance factor A is 
the probability of measuring T > 7^ if the alternative 
hypothesis is true 

/■OO 

A= / p,{T;N)dT. (C2) 

Equivalently, under the frequentist definition of proba- 
bility, it is the fraction of experiments in which the alter- 
native hypothesis is true that measure T > 7^ and thus 
reject the null hypothesis at confidence level R. 



By calculating R and A from the probability distri- 
butions of the test statistic for the null and alternative 
hypothesis as a function of N, an 'acceptance-rejection' 
plot can be built up for each value of and for any given 
level of rejection, Rc, the level of acceptance achiev- 
able for each recoil sample size A'' calculated. In other 
words we find the N such that 'for iV observed recoils, 
the null hypothesis is rejected at the 100iic% confidence 
level in 100Ac% of experiments in which the alternative 
hypothesis is true'. 

Clearly, a high value of Rc is required to reject the 
null hypothesis at high confidence. A high acceptance is 
also required; if, for instance, Ac = 0.1 then only 1 in 10 
experiments will be able to reject the null hypothesis at 
the given Rc, furthermore if Ac is low, the null hypothesis 
could sometimes be erroneously rejected at a given confi- 
dence level with a low number of events due to statistical 
fluctuations. We therefore use Ac = Rc = 0.9 (0.95) as 
OTir criteria, and calculate the corresponding minimum 
number of events required. 
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